
pacman::p_load(tidyverse,here, plm, AER, stargazer,lfe,gtools, permute,parallel,
               doParallel,foreach,cowplot,sp,spdep, splm, igraph,
               did, PanelMatch, ggpubr, fixest, tipr)

#get directory
this_dir <- dirname(rstudioapi::getActiveDocumentContext()$path)

NoLow <- readRDS(file = here(dirname(this_dir), "replication", "data", "main_data_dec2022.rds"))

#------------------------------------------------------------------------
# Descriptive statistics
#-------------------------------------------------------------------------

#-----------------------------------------------
# Figure 1: Characterizing Firms and Jobs in the Sample


# left-hand panel

sec_dist <- NoLow %>%
  dplyr::group_by(sector) %>%
  dplyr::summarise(count_sec = length(unique(primary.name)))

sec_dist <- na.omit(sec_dist)



sec_p<-ggplot(sec_dist, aes(x = reorder(sector, count_sec), y =count_sec)) +
  geom_bar(stat = "identity", alpha = .35) +
  coord_flip() +
  labs(y = "# of MCs Hired", x = NULL) +
  theme_bw() +
  scale_y_continuous(breaks = seq(0,40,5))

sec_p

ggsave(plot = sec_p, 
       filename = here(dirname(this_dir),"replication", "images", "HiringBySector.eps"),
       device = cairo_ps,
       height = 6,
       width = 4)

# right-hand panel

desc_df <- NoLow
desc_df$primary.position <- as.character(desc_df$primary.position)

desc_df$primary.position <- ifelse(desc_df$primary.position == "Board Member", "Board of Directors",
                                   ifelse(desc_df$primary.position == "lobbyist", "Lobbyist", 
                                          ifelse(desc_df$primary.position == "Advisory Board", "Board of Advisors", desc_df$primary.position)))

desc_df$primary.position <- ifelse(desc_df$primary.position == "Consultant" | desc_df$primary.position == "Advisor" | desc_df$primary.position == "Senior Advisor" | desc_df$primary.position == "Outside Consultant" | desc_df$primary.position == "Fund Advisor", "Political Advisor", desc_df$primary.position)

desc_df$primary.position <- ifelse(desc_df$primary.position == "Vice-president" | desc_df$primary.position == "Vice-President" |desc_df$primary.position == "CEO" | desc_df$primary.position == "Secretary" | desc_df$primary.position == "Managin Director", "Management", desc_df$primary.position)

desc_df$primary.position <- ifelse(desc_df$primary.position == "Chairman" | desc_df$primary.position == "Independent Director" | desc_df$primary.position == "Chairman of the Board of Directors", "Board of Directors", desc_df$primary.position)

desc_df$primary.position <- ifelse(desc_df$primary.position == "Director of Federal Affairs" |desc_df$primary.position == "Chief Strategic Officer" | desc_df$primary.position == "Vice Chairman/Sr Policy Advisor", "Political Advisor", desc_df$primary.position)

job_dist <- desc_df %>%
  dplyr::group_by(primary.position) %>%
  dplyr:: summarise(count_job = length(unique(primary.name)))

job_dist <- na.omit(job_dist)


job_p<- ggplot(job_dist, aes(x = reorder(primary.position, count_job), y =count_job)) +
  geom_bar(stat = "identity", alpha = .35, width = .5) +
  coord_flip() +
  labs(y = "# of Former MCs", x = NULL) +
  theme_bw() +
  scale_y_continuous(breaks = seq(0,60,10))


job_p



ggsave(plot = job_p, 
       filename = here(dirname(this_dir),"replication", "images", "HiringByJobType.eps"),
       device = cairo_ps,
       height = 6,
       width = 4)

#--------------------------------------------------------------------
# Results
#------------------------------------------------------------------

# Figure 2: Corporate Tax Rate and Time Until Revolving Door Hire.

time.to.hire <- subset(NoLow, is.na(lead_tax) == F & is.na(tbe) == F)
time.to.hire <- pdata.frame(time.to.hire, index = c("id.num", "year"))

time.to.hire$lead_tax_pct <-  plm::lag(time.to.hire$tax.pct, -1)


low.Lobby <- lowess(x = time.to.hire$tbe, y = time.to.hire$lead_tax_pct, f = 1/4, iter = 3)
low.Lobby2 <- do.call("cbind", low.Lobby)
low.Lobby2 <- as.data.frame(low.Lobby2)


low.plot<- ggplot(time.to.hire, aes(x = -tbe, y = lead_tax_pct)) +
  geom_point() +
  geom_line(data = low.Lobby2, aes(x = -x, y = y)) +
  theme_bw() + theme(panel.grid.major = element_blank(),
                     panel.grid.minor = element_blank(), axis.line = element_line(colour = "black")) +
  scale_x_continuous(limits = c(-5,0),
                     breaks = seq(-5,0), 
                     labels = c("5 Years\nRemaining", "4", "3", "2", "1",
                                "Legislator\nHired")) +
  labs(y = "Effective Tax Rate \n at time t + 1", 
       x = "Years Remaining to Legislator is Hired") +
  scale_y_continuous(limits = c(0,40),
                     breaks = seq(0,40,10),
                     labels = as.character(paste(seq(0,40,10),"%",sep="")))

low.plot

ggsave(plot = low.plot, 
       filename = here(dirname(this_dir),"replication", "images", "lowessPlot.eps"),
       device = cairo_ps,
       height = 4,
       width = 6)

#------------------------------------
# Table 1: The Revolving Door and Corporate Tax Rates


CTRL <- how(within = permute::Within(type = "series"), # time-series ordering ...
            #plots = Plots(strata = NoLow$id.num)
            blocks = NoLow$id.num) # ... within firm-clusters

#create 3000 firm-cluster-permuted treatment assignments

set.seed(185)
perm <- shuffleSet(NoLow$treat2, nset = 3000, control = CTRL)#each row of 'perm' represents a unique assignment


forms_to_est <- list("lead_tax ~ treat2 | id.num + year",
                     "lead_tax ~ treat2 | id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                     "lead_tax ~ treat2 + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year",
                     "lead_tax ~ treat2 + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec)",
                     "lead_tax ~ treat2 + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec) + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)"
)

ri_to_est<- list("lead_tax ~ new_treat | id.num + year",
                 "lead_tax ~ new_treat | id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_tax ~ new_treat + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year",
                 "lead_tax ~ new_treat + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec)",
                 "lead_tax ~ new_treat + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec) + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)"
)

main_est <- rep(NA, length(forms_to_est))
rand_se <- rep(NA, length(forms_to_est))
exact_p <- rep(NA, length(forms_to_est))
list_of_mods <- list()

for(f in 1:length(forms_to_est)){
  
  list_of_mods[[f]] <- felm(formula = as.formula(forms_to_est[[f]]),
                            data = NoLow)
  
  main_est[f] <- coef(list_of_mods[[f]])[1]
  
  
  cl <- makeCluster(detectCores()-1)
  registerDoParallel(cl) # register the cluster
  
  res <- foreach(i = 1:nrow(perm), 
                 .combine = "rbind", 
                 .packages = "lfe") %dopar% {
                   
                   crnt_perm <- perm[i,]
                   NoLow$new_treat <- NoLow$treat2[crnt_perm]
                   
                   this_mod <-felm(formula = as.formula(ri_to_est[[f]])
                                   , data = NoLow)
                   
                   # return the coefficients
                   coef(this_mod)[1]
                 }
  stopCluster(cl)
  
  rand_se[f] <- sd(res)
  res <- as.data.frame(res)
  #res$treat2 <- as.numeric(as.character(res$treat2))
  exact_p[f] <- nrow(as.data.frame(subset(res$new_treat, abs(res$new_treat) > abs(main_est[f]))))/nrow(res)
  
}


NoLow <- pdata.frame(NoLow, index = c("id.num", "year"))


set.seed(185)
plac_perm <- shuffleSet(NoLow$treat_lead, nset = 3000, control = CTRL)#each row of 'perm' represents a unique assignment


plac_to_est <- list("lead_tax ~ treat_lead | id.num + year",
                    "lead_tax ~ treat_lead | id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_tax ~ treat_lead + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year",
                    "lead_tax ~ treat_lead + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec)",
                    "lead_tax ~ treat_lead + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec) + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)"
)

ri_plac_to_est<- list("lead_tax ~ new_plac | id.num + year",
                      "lead_tax ~ new_plac | id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                      "lead_tax ~ new_plac + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year",
                      "lead_tax ~ new_plac + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec)",
                      "lead_tax ~ new_plac + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec) + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)"
)



plac_est <- rep(NA, length(plac_to_est))
rand_se_plac <- rep(NA, length(plac_to_est))
exact_p_plac <- rep(NA, length(plac_to_est))
list_of_plac <- list()

for(f in 1:length(plac_to_est)){
  
  list_of_plac[[f]] <- felm(formula = as.formula(plac_to_est[[f]]),
                            data = NoLow)
  
  plac_est[f] <- coef(list_of_plac[[f]])[1]
  
  
  cl <- makeCluster(detectCores()-1)
  registerDoParallel(cl) # register the cluster
  
  res <- foreach(i = 1:nrow(plac_perm), 
                 .combine = "rbind", 
                 .packages = "lfe") %dopar% {
                   
                   crnt_perm <- plac_perm[i,]
                   NoLow$new_plac <- NoLow$treat_lead[crnt_perm]
                   
                   this_mod <-felm(formula = as.formula(ri_plac_to_est[[f]])
                                   , data = NoLow)
                   
                   # return the coefficients
                   coef(this_mod)[1]
                 }
  stopCluster(cl)
  
  rand_se_plac[f] <- sd(res)
  res <- as.data.frame(res)
  #res$treat2 <- as.numeric(as.character(res$treat2))
  exact_p_plac[f] <- nrow(as.data.frame(subset(res$new_plac, abs(res$new_plac) > abs(plac_est[f]))))/nrow(res)
  
}


# print out table -- make some aesthetic changes by hand in latex.
stargazer(list_of_mods, keep = "treat", 
          se = list(rand_se[1], rand_se[2], rand_se[3], rand_se[4], rand_se[5]),
          covariate.labels = "Revolving Door", omit.stat = c("rsq","adj.rsq",  "f"), 
          df =F, dep.var.labels = "ln Tax Rate",
          column.labels = c("Diff-in-Diff", "MC History", "Firm Covariates", "Firm Sector", "All Adjustments"),
          add.lines = list(c("Exact RI P-value", round(exact_p, digits = 3)),
                           c("Pre-trend", round(plac_est, digits = 3)),
                           c(" ", round(rand_se_plac, digits=3)),
                           c("Firm-level Covariates?", "No", "No", "Yes", "Yes", "Yes"), 
                           c("Firm FE?", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE X MC History", "No", "Yes", "No", "No", "Yes"),
                           c("Year FE X NAICS 2-Digit?", "No", "No", "No", "Yes", "Yes")),
          title = c("The Revolving Door and Corporate Tax Rates"),
          label = c("tab:res"),
          star.cutoffs = NA,
          out = here(dirname(this_dir),"replication", "tables", "table1.tex"))


#------------------------------------------------
# Table 2: audit results

forms_to_est <- list("lead_audit ~ treat2 | id.num + year",
                     "lead_audit ~ treat2 | id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                     "lead_audit ~ treat2 + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year",
                     "lead_audit ~ treat2 + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec)",
                     "lead_audit ~ treat2 + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec) + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)"
)

ri_to_est<- list("lead_audit ~ new_treat | id.num + year",
                 "lead_audit ~ new_treat | id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_audit ~ new_treat + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year",
                 "lead_audit ~ new_treat + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec)",
                 "lead_audit ~ new_treat + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec) + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)"
)

audit_est <- rep(NA, length(forms_to_est))
rand_se_audit <- rep(NA, length(forms_to_est))
exact_p_audit <- rep(NA, length(forms_to_est))
list_of_audit <- list()


for(f in 1:length(forms_to_est)){
  
  list_of_audit[[f]] <- felm(formula = as.formula(forms_to_est[[f]]),
                             data = NoLow)
  
  audit_est[f] <- coef(list_of_audit[[f]])[1]
  
  
  cl <- makeCluster(detectCores()-1)
  registerDoParallel(cl) # register the cluster
  
  res <- foreach(i = 1:nrow(perm), 
                 .combine = "rbind", 
                 .packages = "lfe") %dopar% {
                   
                   crnt_perm <- perm[i,]
                   NoLow$new_treat <- NoLow$treat2[crnt_perm]
                   
                   this_mod <-felm(formula = as.formula(ri_to_est[[f]])
                                   , data = NoLow)
                   
                   # return the coefficients
                   coef(this_mod)[1]
                 }
  stopCluster(cl)
  
  rand_se_audit[f] <- sd(res)
  res <- as.data.frame(res)
  exact_p_audit[f] <- nrow(as.data.frame(subset(res$new_treat, abs(res$new_treat) > abs(audit_est[f]))))/nrow(res)
  
}

set.seed(185)
plac_perm <- shuffleSet(NoLow$treat_lead, nset = 3000, control = CTRL)#each row of 'perm' represents a unique assignment


plac_to_est <- list("lead_audit ~ treat_lead | id.num + year",
                    "lead_audit ~ treat_lead | id.num + year + as.factor(year):as.factor(house) + as.factor(year):as.factor(senate) + as.factor(year):as.factor(governor) + as.factor(year):as.factor(sec) + as.factor(year):factor(n_positions)",
                    "lead_audit ~ treat_lead + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year",
                    "lead_audit ~ treat_lead + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec)",
                    "lead_audit ~ treat_lead + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec)+ as.factor(year):as.factor(house) + as.factor(year):as.factor(senate) + as.factor(year):as.factor(governor) + as.factor(year):as.factor(sec) + as.factor(year):factor(n_positions)"
)

ri_plac_to_est<- list("lead_audit ~ new_plac | id.num + year",
                      "lead_audit ~ new_plac | id.num + year + as.factor(year):as.factor(house) + as.factor(year):as.factor(senate) + as.factor(year):as.factor(governor) + as.factor(year):as.factor(sec) + as.factor(year):factor(n_positions)",
                      "lead_audit ~ new_plac + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year",
                      "lead_audit ~ new_plac + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec)",
                      "lead_audit ~ new_plac + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + as.factor(year):as.factor(id.sec) + as.factor(year):as.factor(house) + as.factor(year):as.factor(senate) + as.factor(year):as.factor(governor) + as.factor(year):as.factor(sec) + as.factor(year):factor(n_positions)"
)



plac_est <- rep(NA, length(plac_to_est))
rand_se_plac <- rep(NA, length(plac_to_est))
exact_p_plac <- rep(NA, length(plac_to_est))
list_of_plac <- list()
res <- rep(NA, nrow(plac_perm))

for(f in 1:length(plac_to_est)){
  
  list_of_plac[[f]] <- felm(formula = as.formula(plac_to_est[[f]]),
                            data = NoLow)
  
  plac_est[f] <- coef(list_of_plac[[f]])[1]
  
  
  for(i in 1:nrow(plac_perm)){
    
    crnt_perm <- plac_perm[i,]
    NoLow$new_plac <- NoLow$treat_lead[crnt_perm]
    
    
    res[i] <- coef(felm(formula = as.formula(ri_plac_to_est[[f]])
                        , data = NoLow))[1]
    
  }
  
  rand_se_plac[f] <- sd(as.numeric(as.character(res)), na.rm = T)
 
}


# print out table -- make some aesthetic changes by hand in latex.
stargazer(list_of_audit, keep = "treat", 
          se = list(rand_se_audit[1], rand_se_audit[2], rand_se_audit[3],
                    rand_se_audit[4], rand_se_audit[5]),
          covariate.labels = "Revolving Door", omit.stat = c("rsq","adj.rsq",  "f"), 
          df =F, dep.var.labels = "Audit Initiation",
          column.labels = c("Diff-in-Diff", "MC History", "Firm Covariates", "Firm Sector", "All Adjustments"),
          add.lines = list(c("Exact P-value", round(exact_p_audit, digits = 3)),
                           c("Pre-trend", round(plac_est, digits = 3)),
                           c(" ", round(rand_se_plac, digits=3)),
                           c("Firm-level Covariates?", "No", "No", "Yes", "Yes", "Yes"), 
                           c("Firm FE?", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE X MC History", "No", "Yes", "No", "No", "Yes"),
                           c("Year FE X NAICS 2-Digit?", "No", "No", "No", "Yes", "Yes")),
          title = c("The Revolving Door and the Probability of Audit Initiation"),
          label = c("tab:audit"),
          star.cutoffs = NA,
          out = here(dirname(this_dir),"replication", "tables", "table2.tex"))

#----------------------------------------------------------------------------------
# Table 3: Firm Strategy and Consequences of Enforcement Changes



CTRL <- how(within = permute::Within(type = "series"), # time-series ordering ...
            blocks = NoLow$id.num) # ... within firm-clusters

#create 3000 firm-cluster-permuted treatment assignments

set.seed(185)
perm <- shuffleSet(NoLow$treat2, nset = 3000, control = CTRL)#each row of 'perm' represents a unique assignment


cl <- makeCluster(detectCores()-1) # create a cluster with 7 cores
registerDoParallel(cl) # register the cluster

res <- foreach(i = 1:nrow(perm), 
               .combine = "rbind", 
               .packages = "lfe") %dopar% {
                 
                 crnt_perm <- perm[i,]
                 NoLow$new_treat <- NoLow$treat2[crnt_perm]
                 
                 mod1 <- felm(log_lead_penalty ~ new_treat|id.num+year|0|id.num,
                              NoLow)
                 
                 mod2 <- felm(log_lead_limit ~ new_treat|id.num+year|0|id.num,
                              NoLow)

                 mod3 <- felm(log_lead_unrec ~ new_treat|id.num+year|0|id.num,
                              NoLow)
                 
                 # return the coefficients
                 this_res <- data.frame(penalty_res = coef(mod1), 
                                        limit_res1 = coef(mod2),
                                        unrec_res = coef(mod3))
                 return(this_res)
               }
stopCluster(cl) # shut down the cluster



mod1 <- felm(log_lead_penalty ~ treat2|id.num+year|0|id.num,
             NoLow)

mod2 <- felm(log_lead_limit ~ treat2|id.num+year|0|id.num,
             NoLow)


mod3 <- felm(log_lead_unrec ~ treat2|id.num+year|0|id.num,
             NoLow)

p1 <- nrow(subset(res, abs(penalty_res) > abs(coef(mod1))) )/3000
p2 <- nrow(subset(res, abs(limit_res1) > abs(coef(mod2)) ))/3000
p3 <- nrow(subset(res, abs(unrec_res) > abs(coef(mod3)) ))/3000

ps <- c(round(p3, digits =2),
        round(p1, digits = 2),
        round(p2, digits = 2) )

stargazer(mod3, mod1, mod2, keep = "treat", 
          se = list(sd(res$penalty_res), sd(res$limit_res1),
                    sd(res$unrec_res)),
          covariate.labels = "Revolving Door", omit.stat = c("rsq","adj.rsq",  "f"), 
          df =F, dep.var.labels = c("ln Unrecognized Benefifts",
                                    "ln Penalties", 
                                    "ln Statue of Limits."),
          column.labels = c(),
          add.lines = list(c("Exact RI P-value", ps),
                           c("Firm FE?", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes")),
          title = c("Consequences of Enforcement Changes"),
          label = c("tab:tax_mech"),
          star.cutoffs = NA,
          out = here(dirname(this_dir),"replication", "tables", "table3.tex"))






#----------------------------------------------------------------------------
# Figure 3: Heterogeneous Effects for Well-Connected Legislators

# functions to compute marginal effects from interactions
marginsplotdf<-function(model, xterm, zterm, zseq){
  coefs<-coef(model)
  cov<-vcovBK(model)
  intterm<-ifelse(is.na(coefs[paste(xterm,zterm,sep=":")]),paste(zterm,xterm,sep=":"),paste(xterm,zterm,sep=":"))
  dy.dx<-coefs[xterm]+coefs[intterm]*zseq
  se.dy.dx<-sqrt(cov[xterm,xterm]+zseq^2*cov[intterm,intterm]+zseq*2*cov[xterm,intterm])
  CI <- data.frame(lwr95 = dy.dx-(1.96*se.dy.dx),
                   lwr90 = dy.dx-(1.68*se.dy.dx),
                   upr90 = dy.dx+(1.68*se.dy.dx),
                   upr95 = dy.dx+(1.96*se.dy.dx))
  
  margins<-data.frame(z=zseq,dydx=dy.dx,se=se.dy.dx, CI)
  return(margins)
}

threeway<-function(model, base, xterm, zterm, xseq, zseq,intterm){
  coefs<-coef(model)
  cov<-vcovBK(model)
  dy.dx <- coefs[base] + # effect of treatment with both moderators at zero
    coefs[xterm]*xseq + # plus additional effect from increasing first moderator ... 
    (coefs[intterm]*zseq) # ... with third moderator at held at informative values
  
  se.dy.dx <- sqrt(cov[xterm,xterm] + 
                     zseq^2*cov[intterm,intterm] + 
                     zseq*2*cov[xterm,intterm])
  
  CI <- data.frame(lwr95 = dy.dx-(1.96*se.dy.dx),
                   lwr90 = dy.dx-(1.68*se.dy.dx),
                   upr90 = dy.dx+(1.68*se.dy.dx),
                   upr95 = dy.dx+(1.96*se.dy.dx))
  
  margins<-data.frame(z = xseq,
                      x = xseq,
                      dydx = dy.dx, 
                      se=se.dy.dx, 
                      CI)
  return(margins)
}

inter_plot <- function(df, z, dydx, x_axis_lab, y_axis_lab, plot_title,
                       x_coord, y_coord, effect_lab,x_breaks,x_labs,
                       ylims, xlims,
                       df2, moderator){
  
  p<-ggplot(df, aes(x = z, y = dydx)) +
    scale_y_continuous(limits = ylims) +
    scale_x_continuous(limits = xlims,
                       breaks = x_breaks,
                       labels = x_labs) +
    #geom_point() +
    geom_line() +
    theme_classic() +
    geom_ribbon(aes(ymin = lwr90, ymax = upr90), alpha = .4) +
    geom_ribbon(aes(ymin = lwr95, ymax = upr95), alpha = .3) +
    geom_hline(yintercept = 0, lty = 2)+
    #geom_hline(yintercept = 0.025, lty = 3)+
    labs(x = x_axis_lab, 
         y =y_axis_lab)    +
    #ggtitle(expression(italic(plot_title))) +
    ggtitle(plot_title) +
    annotate(geom ="text",
             x = x_coord, y = y_coord,
             label = effect_lab,
             size= 3.5)
  
  xhist <- axis_canvas(p, axis = "x") +
    geom_histogram(data = as.data.frame(df2), aes(x = get(moderator)), 
                   fill = "grey") 
  # 
  combined_plot <- insert_xaxis_grob(p, xhist, position = "bottom")
  
}

# calculate models and marginal effects


# average betweeness during tenure
btwn_mod <- plm(lead_tax ~ treat2*av_btwn_con ,
            model = "within", effect = "twoways",
            data = NoLow)

marg_btwn <- marginsplotdf(model = btwn_mod, 
                           xterm = "treat2", 
                           zterm = "av_btwn_con", 
                           zseq = seq(min(NoLow$av_btwn_con, na.rm = T), 
                                      max(NoLow$av_btwn_con, na.rm = T), 
                                      by =.5))


#member of committees with oversight of IRS
over_mod <- plm(lead_tax ~ treat2*TaxOversight ,
            model = "within", effect = "twoways",
            data = NoLow)

marg_over <- marginsplotdf(model = over_mod, 
                           xterm = "treat2", 
                           zterm = "TaxOversight", 
                           zseq = 0:1)


# both combined in three-way interaction
three_mod <- plm(lead_tax ~ treat2*TaxOversight*av_btwn_con,
             model = "within", effect = "twoways",
             data = NoLow)

marg_three <- threeway(model = three_mod, 
                       base = "treat2",
                       intterm = "treat2:TaxOversight:av_btwn_con",
                       xterm = "treat2:av_btwn_con", 
                       zterm = "treat2:TaxOversight", 
                       xseq = seq(min(NoLow$av_btwn_con, na.rm = T), 
                                  max(NoLow$av_btwn_con, na.rm = T), 
                                  by =.5),
                       zseq = 0)

marg_three2 <- marginsplotdf(model = three_mod, 
                             xterm = "treat2", 
                             zterm = "TaxOversight:av_btwn_con", 
                             zseq = seq(min(NoLow$av_btwn_con, na.rm = T), 
                                        max(NoLow$av_btwn_con, na.rm = T), 
                                        by =.5))


cl <- makeCluster(detectCores()-1) # create a cluster with 7 cores
registerDoParallel(cl) # register the cluster

ri_btwn <- foreach(i = 1:nrow(perm), 
                   .combine = "rbind", 
                   .packages = "lfe") %dopar% {
                     
                     crnt_perm <- perm[i,]
                     NoLow$new_treat <- NoLow$treat2[crnt_perm]
                     
                     this_mod <-felm(lead_tax ~ new_treat*av_btwn_con
                                     |id.num + year, data = NoLow)
                     
                     #add the following to restrict sample to av_btwn_con < 2
                     #data = subset(NoLow, av_btwn_con < 2)
                     
                     # return the coefficients
                     coef(this_mod)[3]
                   }
stopCluster(cl)


cl <- makeCluster(detectCores()-1) # create a cluster with 7 cores
registerDoParallel(cl) # register the cluster

ri_oversight <- foreach(i = 1:nrow(perm), 
                        .combine = "rbind", 
                        .packages = "lfe") %dopar% {
                          
                          crnt_perm <- perm[i,]
                          NoLow$new_treat <- NoLow$treat2[crnt_perm]
                          
                          this_mod <-felm(lead_tax ~ new_treat*TaxOversight
                                          |id.num + year, data = NoLow)
                          
                          # return the coefficients
                          coef(this_mod)[3]
                        }
stopCluster(cl)



cl <- makeCluster(detectCores()-1) # create a cluster with 7 cores
registerDoParallel(cl) # register the cluster

three <- foreach(i = 1:nrow(perm), 
                 .combine = "rbind", 
                 .packages = "lfe") %dopar% {
                   
                   crnt_perm <- perm[i,]
                   NoLow$new_treat <- NoLow$treat2[crnt_perm]
                   
                   this_mod <-felm(lead_tax ~ new_treat*TaxOversight*av_btwn_con
                                   |id.num + year, data = NoLow)
                   
                   # return the coefficients
                   coef(this_mod)[7]
                 }
stopCluster(cl)

# plot marginal effects

btwn_plot <- inter_plot(df = marg_btwn, z = z, dydx = dydx,
              x_breaks = -1:3,
              x_labs = c("1 Std.Dev.\n Below", 
                         "0", "1", "2",
                         "3 Std.Dev.\nAbove"),
              x_axis_lab = "Averages Betweenness \n(Standard deviations vs. Mean in\nMC's Average Congress)", 
              y_axis_lab = "Marginal effect\nof Hiring MC",
              plot_title = expression(italic("B: Connectedness of MC")),
              x_coord = 2.5, y_coord = 0.095,
              effect_lab = paste("Estimate =",
                                 round(coef(btwn_mod)[length(coef(btwn_mod))],2),
                                 "\n SE =",
                                 round(sd(ri_btwn), 2)),
              df2 = NoLow, moderator = "av_btwn_con",
              ylims = c(-1.10, 0.50),
              xlims = c(-1.01,3.1))
ggdraw(btwn_plot)


com_plot <- inter_plot(df = marg_over, z = z, dydx = dydx,
                       x_breaks = 0:1,
                       x_labs = c("Not Member", "Member of IRS\nOversight Committee"),
                       x_axis_lab = "Membership of committee with oversight\n(Finance or Ways and Means)", 
                       y_axis_lab = "Marginal effect\nof Hiring MC",
                       plot_title = expression(italic("A: Committee with oversight of IRS")),
                       x_coord = .75, y_coord = 0.008,
                       effect_lab = paste("Estimate =",
                                          round(coef(over_mod)[length(coef(over_mod))],2),
                                          "\n SE =",
                                          round(sd(ri_oversight), 2)),
                       df2 = NoLow, moderator = "TaxOversight",
                       ylims = c(-0.3,0.1),
                       xlims = c(0,1))
ggdraw(com_plot)

# three way, with oversight = 0

OverLow <- subset(NoLow, TaxOversight == 0)

p3 <- inter_plot(df = marg_three, z = z, dydx = dydx,
                 x_breaks = -1:3,
                 x_labs = c("1 Std.Dev.\n Below", 
                            "0", "1", "2",
                            "3 Std.Dev.\nAbove"),
                 x_axis_lab = "Averages Betweenness \n(Standard deviations vs. Mean in\nMC's Average Congress)", 
                 y_axis_lab = "Marginal effect\nof Hiring MC",
                 plot_title = expression(italic("C: Connectedness - No oversight")),
                 x_coord = 2.5, y_coord = 9,
                 effect_lab = ' ',
                 df2 = OverLow, moderator = "av_btwn_con",
                 ylims = c(-1, 0.55),
                 xlims = c(-1.01,3.1))
ggdraw(p3)

# oversight set to 1

OverLow2 <- subset(NoLow, TaxOversight == 1)
p4 <- inter_plot(df = marg_three2, z = z, dydx = dydx,
                 x_breaks = -1:3,
                 x_labs = c("1 Std.Dev.\n Below", 
                            "0", "1", "2",
                            "3 Std.Dev.\nAbove"),
                 x_axis_lab = "Averages Betweenness \n(Standard deviations vs. Mean in\nMC's Average Congress)", 
                 y_axis_lab = "Marginal effect\nof Hiring MC",
                 plot_title = expression(italic("D: Connectedness and oversight")),
                 x_coord = 2.5, y_coord = 0.01,
                 effect_lab = paste("Estimate =",
                                    round(coef(three_mod)[length(coef(three_mod))],2),
                                    "\n SE =",
                                    round(sd(three), 2)),
                 df2 = OverLow2, moderator = "av_btwn_con",
                 ylims = c(-2.5, 0.7),
                 xlims = c(-1.01,3.1))
ggdraw(p4)


top_row <- plot_grid(com_plot, btwn_plot,  nrow = 1)

top_title <- ggdraw() + draw_label("A and B: Twoway interactions", fontface='bold')
top_row <- plot_grid(top_title, top_row, ncol=1, rel_heights=c(0.1, 1)) 


btm_row <- plot_grid(p3, p4, nrow = 1)

btm_title <- ggdraw() + draw_label("C and D: Threeway interaction", fontface='bold')
btm_row <- plot_grid(btm_title, btm_row, ncol=1, rel_heights=c(0.1, 1)) 


final_valid <- plot_grid(top_row, btm_row, ncol =1)
ggdraw(final_valid)


ggsave(final_valid, 
       filename = here(dirname(this_dir),"replication", "images", "InteractionPlots.eps"),
       device = cairo_ps,
       height = 10,
       width = 10)

#---------------------------------------------------------------------------------------
# Figure 4: Hiring a former MC and Other Influence-Seeking Strategies.

NoLow$log_expend <- log(NoLow$expend_Zero+1)
NoLow$log_donate <- log(NoLow$tot+1)

other_strat <- c("lobby~ treat2|id.num+year|0|id.num", 
                 "lobby_irs ~ treat2|id.num+year|0|id.num", 
                 "donate ~ treat2|id.num+year|0|id.num",
                 "other_connex ~ treat2|id.num+year|0|id.num",
                 "irs ~ treat2|id.num+year|0|id.num",
                 "cong~ treat2|id.num+year|0|id.num",
                 "revolver~ treat2|id.num+year|0|id.num")

other_strat_func <- function(DV){
  
  mod <- felm(formula = as.formula(DV),
              data = NoLow,
              Nboot = 500,
              nostat=structure(FALSE, boot= T),
              bootcluster = factor(NoLow$id.num))
  
  reps <- mod$boot
  reps <- as.data.frame(reps)
  #CI <- quantile(reps$reps, probs =c(0.025,0.05,0.95,0.975))
  
  
  res <- data.frame(est = coef(mod), se = sd(reps$reps))
  res <- res %>%
    mutate(lwr95 = est - 1.96*se,
           lwr90 = est - 1.645*se,
           upr90 = est + 1.645*se,
           upr95 = est + 1.96*se,)
  
  #names(res)<-c("est", "lwr95", "lwr90", "upr90", "upr95")
  return(res)
}


other_strat_res <- map_df(other_strat,
                          ~ other_strat_func(.))

other_strat_res$DV <- c("Lobbying", "Lobbying the IRS", 
                        "Campaign Donations", "Hire Staffer/IRS/CEA",
                        "Hire IRS", "Hire Staffer", "Revolver on Contract")

p <- ggplot(other_strat_res, aes(x = est, y = DV)) +
  geom_point() +
  geom_errorbarh(aes(xmin = lwr95,
                     xmax = upr95), height = 0, colour = "grey") +
  geom_errorbarh(aes(xmin = lwr90,
                     xmax = upr90), height = 0, colour = "black") +
  theme_classic() +
  geom_vline(xintercept = 0, lty = 2) +
  labs(x = "Correlation between Hiring MC and\nOther Strategies",
       y = NULL) +
  scale_x_continuous(limits = c(-0.1, 0.1))

p

ggsave(filename = here(dirname(this_dir),"replication", "images", "OtherStrategies.eps"),
       device = cairo_ps,
       width = 6.5, height = 6.5)



#-------------------------------------------------------------------
# APPENDIX
#------------------------------------------------------------------------

#-----------------------------------------
# Appendix B: Descriptive statistics

# Table B.1: Firm Finances

desc <- select(NoLow, treat2, tax.pct, total_capital, DS_total_assets, enterprise_value, employees,
               net_sales_revenue, gross_income, turnover_volume, market_value, prices_adj,
)
desc <- as.data.frame(desc)

stargazer(desc, summary.stat = c("n", "median", "sd", "p25","p75"),
          covariate.labels = c("Years w/ Revolver (Proportion)", 
                               "Effective Tax Rate (%)",
                               "$ Total Capital (Thousands)",
                               "$ Total Assets (Thousands)",
                               "$ Enterprise Value (Thousands)",
                               "# Employees (Thousands)",
                               "$ Net Sales (Thousands)",
                               "$ Gross Income (Thousands)",
                               "Turnover Volume (# Shares)",
                               "$ Markel Value (Millions)",
                               "$ Stock Price"),
          dep.var.labels = c("Mean", "Std. Dev.", "25th Pctile", "75th Pctile"),
          out = here(dirname(this_dir),"replication", "tables", "firm_descriptives.tex"))

# the table includes the mean of the treatment variable and the median for the 
# rest of the variables.
# the following produces that, which can then be inputted into the table manually
mean(NoLow$treat2, na.rm=T)


#Table B.2: Legislator Characteristics and Firm Political Behavior

mc_desc <- select(NoLow, av_btwn_con, TaxOversight, agency_expend)
mc_desc <- as.data.frame(mc_desc)
mc_desc$agency_expend <- ifelse(is.na(mc_desc$agency_expend)==T, 0,
                                mc_desc$agency_expend)

stargazer(mc_desc, summary.stat = c("n", "mean", "sd", "p25","p75"),
          covariate.labels = c("Average Betweenness", "Membership of Tax Committee",
                               "Lobbying the IRS"),
          dep.var.labels = c("Mean", "Std. Dev.", "25th Pctile", "75th Pctile"),
          out = here(dirname(this_dir),"replication", "tables", "legislator_descriptives.tex"))

#Figure B.1: How Many Directorships Do former Legislators Hold?

sub_desc <- subset(desc_df, primary.position == "Board of Directors")


dir_dist <- sub_desc %>%
  dplyr::group_by(primary.name) %>%
  dplyr::summarise(count_dir = length(unique(id.num)))


dir_p <- ggplot(dir_dist, aes(x = count_dir)) +
  geom_density(fill = "black", alpha = .2, colour = "white") +
  labs(y = "Density", x = "# of Directorships Held by Former MC") +
  theme_bw() +
  scale_x_continuous(breaks = 1:7)

dir_p

ggsave(plot = dir_p,
       filename = here(dirname(this_dir),"replication", "images", "DistDirectorships.eps"),
       device = cairo_ps,
       width = 6.5, height = 5)

#-----------------------------------------------------
# Appendix C: Extracting Information on IRS Audits from 10-Ks

#Figure C.1: How Often are Firms Audited by the IRS?

aud_dist <- desc_df %>%
  dplyr::group_by(id.num) %>%
  dplyr::summarise(n_audit = sum(IRS_audit, na.rm = T))

aud_p <- ggplot(aud_dist, aes(x = n_audit)) +
  geom_density(fill = "black", alpha = .2, colour = "white") +
  labs(y = "Density", x = "# of Audits Throughout Period") +
  theme_bw() +
  scale_x_continuous(breaks = seq(0,10,2))

aud_p
ggsave(plot = aud_p,
       filename = here(dirname(this_dir),"replication", "images", "DistAudits.eps"),
       device = cairo_ps,
       width = 6.5, height = 5)

#------------------------------------------------------------------------
#D Additional Analyses

#-----------------------
# D.1 Correlation with Firm Economic Characteristics

#Table D.1: No Correlation Between Hiring and Changes in Firm Characteristics
full.bal <- plm(treat2 ~ lag_tax +
                  lag_capital + 
                  lag_assets +
                  lag_ent +
                  lag_emp +
                  lag_sales +
                  lag_income +
                  lag_turnover +  
                  lag_market_value +  
                  lag_prices,
                  #factor(year), 
                data = NoLow, model = "within", effect = "twoways")

se.full.bal <- sqrt(diag(vcovBK(full.bal, type = "HC1", cluster = "group")))
coeftest(full.bal, vcov = vcovBK(full.bal, type = "HC1", cluster = "group"))

w_test <- pwaldtest(full.bal, vcov=vcovHC, type = "HC1")
as.vector(round(w_test$statistic, digits = 3))
as.vector(round(w_test$p.value, digits = 3))


stargazer(full.bal,
          se = list(se.full.bal),
          covariate.labels = c("Tax Rate$_{t-1}$", 
                               "ln Total Capital$_{t-1}$",
                               "ln Total Assets$_{t-1}$", 
                               "ln Enterprise Value$_{t-1}$", 
                               "ln Employees$_{t-1}$", 
                               "ln Net Revenue$_{t-1}$",
                               "ln Gross Income$_{t-1}$", 
                               "ln Turnover Volume$_{t-1}$", 
                               "ln Market Value$_{t-1}$", 
                               "Share Price$_{t-1}$"),
          add.lines = list(c("Wald Stat", as.vector(round(w_test$statistic, digits = 3))),
                           c("Wald P value", as.vector(round(w_test$p.value, digits = 3))),
                           c("Firm FEs?", "Yes"),
                           c("Time FEs?", "Yes")),
          omit.stat = "f", no.space = T, font.size = "footnotesize",
          out = here(dirname(this_dir),"replication", "tables", "balance.tex"))


#----------------------
# D.2 Placebo Sample of Foreign Firms

# Table D.2: Results for Sample of Foreign Firmsn

non_us <- subset(NoLow, ccy != "U$")

#identify clusters
CTRL <- how(within = permute::Within(type = "series"), # time-series ordering ...
            plots = Plots(strata = non_us$id.num)) # ... within firm-clusters

set.seed(185)
us_perm <- shuffleSet(non_us$treat2, nset = 3000, control = CTRL)#each row of 'perm' represents a unique assignment


cl <- makeCluster(detectCores()-1) # create a cluster with 7 cores
registerDoParallel(cl) # register the cluster

res <- foreach(i = 1:nrow(us_perm), 
               .combine = "rbind", 
               .packages = "lfe") %dopar% {
                 
                 crnt_perm <- us_perm[i,]
                 non_us$new_treat <- non_us$treat2[crnt_perm]
                 
                 this_mod <-felm(lead_tax ~ new_treat 
                                 |id.num + year, data = non_us)
                 
                 # return the coefficients
                 coef(this_mod)[1]
               }
stopCluster(cl) # shut down the cluster


res <- as.data.frame(res)

us_mod <- felm(lead_tax ~ treat2 
               | id.num + year ,
               data = non_us) 

us_est <-coef(us_mod)[1]

crit <- nrow(as.data.frame(subset(res$new_treat, abs(res$new_treat) > abs(us_est))))/nrow(res)

se_us <- sd(res$new_treat)
stargazer(us_mod, covariate.labels = "Revolving Door", 
          se = list(se_us),
          dep.var.labels = "Tax Rate", df=F,
          omit.stat = c("rsq","adj.rsq",  "f"),
          add.lines = list(c("Exact P-Value", round(crit, digits = 3)),
                           c("Firm FE?", "Yes"),
                           c("Year FE?", "Yes")),
          out = here(dirname(this_dir),"replication", "tables", "foreign_firms.tex"))

#--------------------------------------------------------------------------
# E Robustness Checks

#-----------------------------
# E.1 Controlling for Pre-Hiring Tax Rates


NoLow <- pdata.frame(NoLow, index = c("id.num", "year"))


biv.model <- plm(lead_tax ~ lag_tax +
                   treat2, 
                 data = NoLow, model = "within", effect = "twoways")
se.biv <- sqrt(diag(vcovBK(biv.model, type = "HC1", cluster = "group")))

#controlling for assets
asset.model <- plm(lead_tax ~ lag_tax +
                     treat2 +
                     
                     lag_capital + 
                     lag_assets +
                     lag_ent +
                     lag_emp,
                   data = NoLow, model = "within", effect = "twoways")
se.asset <- sqrt(diag(vcovBK(asset.model, type = "HC1", cluster = "group")))

# controlling for general performance

perf.model <- plm(lead_tax ~ lag_tax +
                    treat2 +
                    
                    lag_capital + 
                    lag_assets +
                    lag_ent +
                    lag_emp+
                    lag_sales +
                    lag_income,
                  data = NoLow, model = "within", effect = "twoways")
se.perf <- sqrt(diag(vcovBK(perf.model, type = "HC1", cluster = "group")))

# add attention on market -- full set of controls

full.model <- plm(lead_tax ~ lag_tax +
                    treat2 +
                    lag_capital + 
                    lag_assets +
                    lag_ent +
                    lag_emp+
                    lag_sales +
                    lag_income +
                    lag_turnover +  
                    lag_market_value +  
                    lag_prices,
                  data = NoLow, model = "within", effect = "twoway")

se.full <- sqrt(diag(vcovBK(full.model, type = "HC1", cluster = "group")))



# placebo models

biv.plac <- plm(plm::lag(log_tax,1) ~ plm::lag(tax.pct,2) +
                  treat2, 
                data = NoLow, model = "within", effect = "twoways")

se.biv.plac <- sqrt(diag(vcovHC(biv.plac, type = "HC0", cluster = "group")))

#controlling for assets
asset.plac <- plm(plm::lag(log_tax,1) ~ plm::lag(tax.pct,2) +
                    treat2 +

                    lag_capital + 
                    lag_assets +
                    lag_ent +
                    lag_emp, 
                  data = NoLow, model = "within", effect = "twoways")
se.asset.plac <- sqrt(diag(vcovHC(asset.plac, type = "HC0", cluster = "group")))

perf.plac <- plm(plm::lag(log_tax,1) ~ plm::lag(tax.pct,2) +
                   treat2 +
                   
                   lag_capital + 
                   lag_assets +
                   lag_ent +
                   lag_emp+
                   lag_sales +
                   lag_income, 
                 data = NoLow, model = "within", effect = "twoways")
se.perf.plac <- sqrt(diag(vcovHC(perf.plac, type = "HC0", cluster = "group")))

plac_full <- plm(plm::lag(log_tax,1) ~ plm::lag(tax.pct,2) +
                   treat2 + 
                   lag_capital + 
                   lag_assets +
                   lag_ent +
                   lag_emp+
                   lag_sales +
                   lag_income +
                   lag_turnover +  
                   lag_market_value +  
                   lag_prices, 
                 data = NoLow, model = "within", effect = "twoways")
se.full.plac <- sqrt(diag(vcovHC(plac_full, type = "HC0", cluster = "group")))


### set up regression table

se_plac <- c(paste("(", round(se.biv.plac[2], digits = 3), ")", sep = ""), 
             paste("(", round(se.asset.plac[2], digits = 3), ")", sep = ""), 
             paste("(", round(se.perf.plac[2], digits = 3), ")", sep = ""), 
             paste("(", round(se.full.plac[2], digits = 3), ")", sep = ""))
pe_plac <- c(round(coef(biv.plac)[2],digits = 3), round(coef(asset.plac)[2], digits = 3),
             round(coef(perf.plac)[2], digits = 3), round(coef(plac_full)[2], digits = 3))

stargazer(biv.model, asset.model, perf.model,full.model,
          #biv.plac, asset.plac, perf.plac, plac_full,
          covariate.labels = c("Tax Rate$_{t-1}$", 
                               "Revolving Door$_{t}$", 
                               "ln Total Capital$_{t-1}$",
                               "ln Total Assets$_{t-1}$", 
                               "ln Enterprise Value$_{t-1}$", 
                               "ln Employees$_{t-1}$", 
                               "ln Net Revenue$_{t-1}$",
                               "ln Gross Income$_{t-1}$", 
                               "ln Turnover Volume$_{t-1}$", 
                               "ln Market Value$_{t-1}$", 
                               "Share Price$_{t-1}$"),
          column.labels = c("Bivariate", "Assets", "Performance", "Attention", "All"),
          se = list(se.biv, se.asset, se.perf, se.full),
          #se.biv.plac, se.asset.plac, se.perf.plac,
          #se.full.plac),
          omit.stat = c("adj.rsq", "f", "rsq"), no.space = T, font.size = "footnotesize",
          header= F,  
          dep.var.labels = c("ln Effective Tax Rate$_{t+1}$"), 
          df = F,
          out = here(dirname(this_dir),"replication", "tables", "control_for_lag_DV.tex"),
          title = "The Revolving Door and Corporate Tax Rates",
          add.lines = list(c("Pre-Trend", pe_plac),
                           c(" ", se_plac),
                           c("Company FEs?", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FEs?", "Yes", "Yes", "Yes", "Yes", "Yes")))


#-----------------------------
# E.2 Robustness to Outlying Observations

# Figure E.1: Robustness of the Results to Including Extreme Observations

pan.dat <- readRDS(file = here(dirname(this_dir), "replication", "data", "full_data.rds"))

#pan.dat <- pdata.frame(full_data, index = c("id.num", "year"))

panCTRL <- how(within = permute::Within(type = "series"), # time-series ordering ...
               blocks = pan.dat$id.num) # ... within firm-clusters


set.seed(185)
pan_perm <- shuffleSet(pan.dat$treat2, nset = 3000, control = panCTRL)#each row of 'perm' represents a unique assignment


cl <- makeCluster(detectCores()-1) # create a cluster with 7 cores
registerDoParallel(cl) # register the cluster

res <- foreach(i = 1:nrow(perm), 
               .combine = "rbind", 
               .packages = "lfe") %dopar% {
                 
                 crnt_perm <- pan_perm[i,]
                 pan.dat$new_treat <- pan.dat$treat2[crnt_perm]
                 
                 this_mod <-felm(lead_tax ~ new_treat 
                                 |id.num + year, data = pan.dat)
                 
                 # return the coefficients
                 coef(this_mod)[1]
               }
stopCluster(cl) # shut down the cluster


cl <- makeCluster(detectCores()-1) # create a cluster with 7 cores
registerDoParallel(cl) # register the cluster

base <- foreach(i = 1:nrow(perm), 
                .combine = "rbind", 
                .packages = "lfe") %dopar% {
                  
                  crnt_perm <- perm[i,]
                  NoLow$new_treat <- NoLow$treat2[crnt_perm]
                  
                  this_mod <-felm(lead_tax ~ new_treat 
                                  |id.num + year, data = NoLow)
                  
                  # return the coefficients
                  coef(this_mod)[1]
                }
stopCluster(cl) 

sample_est <- coef( felm(lead_tax ~ treat2 
                         | id.num + year ,
                         data = pan.dat) )[1]

baseline_est <- coef( felm(lead_tax ~ treat2| id.num + year, data = NoLow) )[1]

est <- data.frame(pe = rbind(sample_est, baseline_est), 
                  se = rbind(sd(res), sd(base)))

est$spec <- c("All Observations", "Baseline\n(Excl. top/bottom 2.5%)")

excl <- ggplot(est, aes(x = treat2, y = spec)) +
  geom_point() +
  geom_errorbarh(aes(xmin = treat2 - (1.96*se),
                     xmax = treat2 + (1.96*se)), height=0) +
  theme_classic() +
  labs(x = "Coefficient on Revolving Door", y = NULL) +
  geom_vline(xintercept = 0, lty = 3)

excl

ggsave(plot = excl,
       filename = here(dirname(this_dir),"replication", "images", "RobustToExcl.eps"),
       device = cairo_ps,
       width = 5, height = 6.5)



# Figure E.2: Distribution of Leave-One-Out Estimates

jack <- rep(NA, nrow(NoLow))

for(i in 1:nrow(NoLow)){
  
  sub_low <- NoLow[-i, ]
  
  first <- plm(plm::lag(log_tax,-1) ~ #lag_tax +
                 treat2,
               data = sub_low, model = "within", effect = "twoways")
  jack[i] <- coef(first)[1]
  
}

jack <- data.frame(jack, NoLow$id.num)

jack_dist<-ggplot(jack, aes(x = jack)) +
  geom_histogram(fill = "black", alpha = .25, bins=100) +
  #geom_density() +
  theme_classic() +
  scale_x_continuous(limits = c(-0.1, 0)) +
  geom_vline(xintercept = 0, lty = 3) 

jack_dist

ggsave(jack_dist,
       filename = here(dirname(this_dir),"replication", "images", "JackknoifedDist.eps"),
       device = cairo_ps, width = 6.5, height = 5.1)



#-----------------------------
# E.3 Robustness to Choice of Transformation

OnlyUS <- subset(NoLow, ccy == "U$")

OnlyUS$bi_tax <-  sign(OnlyUS$tax.pct)*log(1+abs(OnlyUS$tax.pct)/10)
OnlyUS$lead_bi_tax <- plm::lag(OnlyUS$bi_tax, -1)

# the inverse hyperbolic sine
OnlyUS$ihs_tax <- log(OnlyUS$tax.pct + (sqrt((OnlyUS$tax.pct^2)+1)))
OnlyUS$lead_ihs_tax <- plm::lag(OnlyUS$ihs_tax, -1)

OnlyUS$lead_level <- plm::lag(OnlyUS$tax.pct,-1)


inter_forms <- c("lead_level ~ treat2| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_level ~ treat2 + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_ihs_tax ~ treat2| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_ihs_tax ~ treat2 + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_bi_tax ~ treat2| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_bi_tax ~ treat2 + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_level ~ treat2*TaxOversight| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_level ~ treat2*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_level ~ treat2*TaxOversight*av_btwn_con| id.num + year +year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_ihs_tax ~ treat2*TaxOversight| id.num + year +year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_ihs_tax ~ treat2*av_btwn_con| id.num + year +year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_ihs_tax ~ treat2*TaxOversight*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_bi_tax ~ treat2*TaxOversight| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_bi_tax ~ treat2*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                 "lead_bi_tax ~ treat2*TaxOversight*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)")

ri_inter_forms <- c("lead_level ~ new_treat| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_level ~ new_treat + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices | id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_ihs_tax ~ new_treat| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_ihs_tax ~ new_treat + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_bi_tax ~ new_treat| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_bi_tax ~ new_treat + lag_capital + lag_assets + lag_ent + lag_emp + lag_sales + lag_income + lag_turnover + lag_market_value + lag_prices| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_level ~ new_treat*TaxOversight| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_level ~ new_treat*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_level ~ new_treat*TaxOversight*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_ihs_tax ~ new_treat*TaxOversight| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_ihs_tax ~ new_treat*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_ihs_tax ~ new_treat*TaxOversight*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_bi_tax ~ new_treat*TaxOversight| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_bi_tax ~ new_treat*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)",
                    "lead_bi_tax ~ new_treat*TaxOversight*av_btwn_con| id.num + year + year:house + year:senate + year:governor + year:sec + year:factor(n_positions)")



#identify clusters
CTRL <- how(within = permute::Within(type = "series"), # time-series ordering ...
            #plots = Plots(strata = NoLow$id.num)
            blocks = OnlyUS$id.num) # ... within firm-clusters

#create 3000 firm-cluster-permuted treatment assignments

set.seed(185)
perm <- shuffleSet(OnlyUS$treat2, nset = 3000, control = CTRL)#each row of 'perm' represents a unique assignment


crnt_mod <- list()
sample_est <- matrix(NA, nrow = length(inter_forms), ncol = 7)
exact_p <-  matrix(NA, nrow = length(inter_forms), ncol = 7)
sd_res <- matrix(NA, nrow = length(inter_forms), ncol = 7)

for(f in 1:6){
  
#for(f in 1:length(inter_forms)){
  
  print(f)
  
  crnt_mod[[f]] <- felm(as.formula(inter_forms[f]),
                        data = OnlyUS)
  crnt_coef <- as.data.frame(coef( crnt_mod[[f]] ) )
  crnt_coef <- crnt_coef[ grep(x= rownames(crnt_coef), pattern = "treat"),]
  
  sample_est[f, 1:length( crnt_coef ) ] <- t(crnt_coef )
  
  cl <- makeCluster(detectCores()-1)
  registerDoParallel(cl) # register the cluster
  
  res <- foreach(i = 1:nrow(perm), 
                 .combine = "rbind", 
                 .packages = "lfe") %dopar% {
                   
                   crnt_perm <- perm[i,]
                   OnlyUS$new_treat <- OnlyUS$treat2[crnt_perm]
                   
                   this_mod <-felm(formula = as.formula(ri_inter_forms[f])
                                   , data = OnlyUS)
                   
                   # return the coefficients
                   coef(this_mod)
                 }
  stopCluster(cl)
  
  
  for(c in 1:length(crnt_coef)){
    
    sd_res[f,c] <- sd(res[,c])
    
    exact_p[f, c] <- nrow(as.data.frame(subset(res[,c], abs(res[,c]) > abs(sample_est[f,c]))))/nrow(res)
    
  }
  
  
}



stargazer(crnt_mod[1:6], keep = "treat", 
          se = list(sd_res[1,], sd_res[2,], sd_res[3,]),
          covariate.labels = "Revolving Door", omit.stat = c("rsq","adj.rsq",  "f"), 
          df =F, 
          dep.var.labels = c("Tax Rate", "Bi-Symlog(Tax Rate)", "IHS(Tax Rate)"),
          column.labels = c("Levels", "Levels", "Bi-Sym", "Bi-Sym", "IHS",  "IHS"),
          add.lines = list(c("Exact P-value", round(exact_p, digits = 3)),
                           c("Covariates?", "No", "Yes", "No", "Yes", "No", "Yes"),
                           c("Firm FE?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year X Chamber FE?", "Yes", "Yes", "Yes", "Yes", "Yes")),
          title = c("Robustness to Choice of Transformation"),
          label = c("tab:RobMain"),
          out = here(dirname(this_dir),"replication", "tables", "Alternative_transform_1.tex"))



crnt_mod <- list()
sample_est <- matrix(NA, nrow = length(inter_forms), ncol = 7)
exact_p <-  matrix(NA, nrow = length(inter_forms), ncol = 7)
sd_res <- matrix(NA, nrow = length(inter_forms), ncol = 7)

inter_forms <- inter_forms[7:15]
ri_inter_forms <- ri_inter_forms[7:15]


for(f in 1:length(inter_forms)){
  
  print(f)
  
  crnt_mod[[f]] <- felm(as.formula(inter_forms[f]),
                        data = OnlyUS)
  crnt_coef <- as.data.frame(coef( crnt_mod[[f]] ) )
  crnt_coef <- crnt_coef[ grep(x= rownames(crnt_coef), pattern = "treat"),]
  
  sample_est[f, 1:length( crnt_coef ) ] <- t(crnt_coef )
  
  cl <- makeCluster(detectCores()-1)
  registerDoParallel(cl) # register the cluster
  
  res <- foreach(i = 1:nrow(perm), 
                 .combine = "rbind", 
                 .packages = "lfe") %dopar% {
                   
                   crnt_perm <- perm[i,]
                   OnlyUS$new_treat <- OnlyUS$treat2[crnt_perm]
                   
                   this_mod <-felm(formula = as.formula(ri_inter_forms[f])
                                   , data = OnlyUS)
                   
                   # return the coefficients
                   coef(this_mod)
                 }
  stopCluster(cl)
  
  
  for(c in 1:length(crnt_coef)){
    
    sd_res[f,c] <- sd(res[,c])
    
    exact_p[f, c] <- nrow(as.data.frame(subset(res[,c], abs(res[,c]) > abs(sample_est[f,c]))))/nrow(res)
    
  }
  
  
}

stargazer(crnt_mod, keep = "treat", 
          se = list(sd_res[1, c(1, NA, 2) ], 
                    sd_res[2, c(1, NA, 2)],
                    sd_res[3, c(1, NA, NA,  2, 3, NA, 4)],
                    sd_res[4, c(1, NA, 2) ], sd_res[5, c(1, NA, 2) ], sd_res[6,c(1, NA, NA,  2, 3, NA, 4)],
                    sd_res[7,c(1, NA, 2)], sd_res[8,c(1, NA, 2)], sd_res[9,c(1, NA, NA,  2, 3, NA, 4)]),
          covariate.labels = c("Revolving Door", "Revolving Door X Oversight",
                               "Revolving Door X Betweenness",
                               "Revolving Door X Oversight X Betweenness"), 
          omit.stat = c("rsq","adj.rsq",  "f"), 
          df =F, dep.var.labels = "Audit Initiation",
          column.labels = c("Levels", "Levels", "Levels", "Bi-Sym", "Bi-Sym", "Bi-Sym", "IHS", "IHS",  "IHS"),
          add.lines = list(c("Firm FE?", "Yes", "Yes", "Yes", "Yes", "Yes","Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes","Yes", "Yes", "Yes", "Yes"),
                           c("Year X Chamber FE?", "Yes", "Yes", "Yes", "Yes", "Yes","Yes", "Yes", "Yes", "Yes")),
          title = c("Robustness of Interactions to Choice of Transformation"),
          label = c("tab:RobInter"),
          out = here(dirname(this_dir),"replication", "tables", "Alternative_transform_2.tex"))



#-----------------------------
# E.4 Intensive Margin of Lobbying and Donations as Alternative Strategy

#Table E.4: Limited Evidence of Relationship between Hiring Revolvers and Intensive Margin of Political Presence

expend_mod <- felm(log(expend_Zero+1) ~ treat2|id.num+year|0|id.num,
                   data = NoLow)

donate_mod <- felm(log(tot+1)~ treat2|id.num+year|0|id.num,
                   data = NoLow)

stargazer(donate_mod, expend_mod,
          covariate.labels = "Hire Revolver",
          dep.var.labels = c("ln Donated Amount",
                             "ln Lobbied Amount"),
          omit.stat = c("rsq", "adj.rsq", "f"),
          df = FALSE,
          title = "Limited Evidence of Relationship between Hiring Revolvers and Intensive Margin of Political Presence",
          out = here(dirname(this_dir),"replication", "tables", "ExtensiveMargin.tex"))


#-----------------------------
# E.5 Alternative Standard Errors

# Table E.5: Results with Abadie et al Clustered Standard Errors

mod <- felm(lead_tax ~ treat2|id.num+year|0|id.num, data = NoLow)
mod2 <- felm(lead_audit ~ treat2|id.num+year|0|id.num, data = NoLow)

# function for estimating the SEs

aaiw_est <- function(df, DV, mod_obj){
  
  #listwise delete
  df <- df %>%
    drop_na({{DV}}, treat2)
  
  # mean under treatment
  correct1 <- df%>%
    filter(treat2 == 1) %>%
    group_by(id.num) %>%
    dplyr::summarise(bar_e_1 = mean({{DV}}),
                     m_1 = n())
  
  #mean under control
  correct0 <- df %>%
    filter(treat2 == 0) %>%
    group_by(id.num) %>%
    dplyr::summarise(bar_e_0 = mean({{DV}}),
                     m_0 = n())
  
  
  # merge together 
  correct <- left_join(correct1, correct0, by = "id.num")
  
  # get total number of observations
  correct$tot <-((correct$m_1 + correct$m_0)) 
  
  # get within firm heterogeneity
  correct$tau_c <- correct$bar_e_1 - correct$bar_e_0
  correct$tau <- coef(mod)[1] # overall different btwn POs
  
  # the correction is a weighted average of the squared difference between firm heterogeneity
  # and the overall average
  correction <- (1/(nrow(df)^2))*sum(correct$tot^2*(correct$tau_c - correct$tau)^2)
  
  # make the correction
  AAIW_se <- sqrt(diag(vcov(mod)) - correction)
  
  return(AAIW_se)
  
}

# get standard errors  
corr_se <- aaiw_est(df=NoLow, DV = lead_tax, mod_obj = mod)
corr_se_audit <- aaiw_est(df=NoLow, DV = lead_audit, mod_obj = mod2)

# t stats
coef(mod)/corr_se
coef(mod2)/corr_se_audit

# upper limits of confidence intervals
coef(mod) + 1.96*corr_se
mean(NoLow$tax.pct)*(coef(mod) + 1.96*corr_se)
mean(NoLow$tax.pct)*(coef(mod))

coef(mod2) + 1.96*corr_se

# how much money is that?
(median(NoLow$pretax_income, na.rm=T)*1000)*(0.2392209/100)
(median(NoLow$pretax_income, na.rm=T)*1000)*(1.365033/100)

# make table

stargazer(mod, mod2, 
          covariate.labels = c("Revolving Door"),
          se = list(corr_se, corr_se_audit), keep.stat ="n",
          add.lines = list(c("Firm FE?", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes")),
          dep.var.labels = c("ln ETR", "Audit"),
          title = "Results with Abadie et al Clustered Standard Errors",
          label = "tab:se",
          out = here(dirname(this_dir),"replication", "tables", "AbadieSE.tex"))




#-----------------------------
# E.6 Missing Data on Covariates

# Table E.6: Does Missingness Correlate with Revolver Arrival?

NoLow$miss_capital <-ifelse(is.na(NoLow$total_capital) == T, 1, 0)
NoLow$miss_assets <- ifelse(is.na(NoLow$DS_total_assets) == T, 1, 0)
NoLow$miss_ent <- ifelse(is.na(NoLow$enterprise_value) == T, 1, 0)
NoLow$miss_emp <- ifelse(is.na(NoLow$employees) == T, 1, 0) 
NoLow$miss_sales <-  ifelse(is.na(NoLow$net_sales_revenue) == T, 1, 0) 
NoLow$miss_income <-  ifelse(is.na(NoLow$gross_income) == T, 1, 0)
NoLow$miss_turnover <- ifelse(is.na(NoLow$turnover_volume) == T, 1, 0)  
NoLow$miss_market_value <- ifelse(is.na(NoLow$market_value) == T, 1, 0)  
NoLow$miss_prices<- ifelse(is.na(NoLow$prices_adj) == T, 1, 0)

miss_dat <- NoLow %>% select(treat2, miss_capital, miss_assets, miss_ent,
                             miss_emp, miss_sales, miss_income, miss_turnover,
                             miss_market_value)

miss_mod <- lm(treat2 ~ ., 
               data = miss_dat)

miss_mod_logit <- glm(treat2 ~ ., 
                      data = miss_dat, family = "binomial")


stargazer(miss_mod, miss_mod_logit,
          covariate.labels = c("Missing Capital", "Missing Assets",
                               "Missing Enterprise Value", "Missing Employees",
                               "Missing Sales", "Missing Income", "Missing Turnover",
                               "Missing Market Value"),
          dep.var.labels = "Revolver", title = "Does Missingness Correlate with Revolver Arrival?",
          keep.stat = c("n", "rsq", "f"), df = FALSE,
          out = here(dirname(this_dir),"replication", "tables", "missing_data.tex"))



#----------------------------------------------------------------------------
# Appendix F Interaction with Legislative Effectiveness Scores


marginsplotdf2<-function(model, xterm, zterm, zseq){
  coefs<-coef(model)
  cov<-vcovHC(model)
  intterm<-ifelse(is.na(coefs[paste(xterm,zterm,sep=":")]),paste(zterm,xterm,sep=":"),paste(xterm,zterm,sep=":"))
  dy.dx<-coefs[xterm]+coefs[intterm]*zseq
  se.dy.dx<-sqrt(cov[xterm,xterm]+zseq^2*cov[intterm,intterm]+zseq*2*cov[xterm,intterm])
  CI <- data.frame(lwr95 = dy.dx-(1.96*se.dy.dx),
                   lwr90 = dy.dx-(1.68*se.dy.dx),
                   upr90 = dy.dx+(1.68*se.dy.dx),
                   upr95 = dy.dx+(1.96*se.dy.dx))
  
  margins<-data.frame(z=zseq,dydx=dy.dx,se=se.dy.dx, CI)
  return(margins)
}


# Figure F.1: Interaction with Legislative Effectiveness Scores. 

CTRL <- how(within = permute::Within(type = "series"), # time-series ordering ...
            blocks = NoLow$id.num) # ... within firm-clusters

set.seed(185)
perm <- shuffleSet(NoLow$treat2, nset = 3000, control = CTRL)#each row of 'perm' represents a unique assignment



les_mod <- plm(lead_tax ~ treat2*av_les, model = "within", effects = "twoway", data = NoLow)

cl <- makeCluster(detectCores()-1)
registerDoParallel(cl) # register the cluster

res <- foreach(i = 1:nrow(perm), 
               .combine = "rbind", 
               .packages = "lfe") %dopar% {
                 
                 crnt_perm <- perm[i,]
                 NoLow$new_treat <- NoLow$treat2[crnt_perm]
                 
                 this_mod <-felm(lead_tax ~ new_treat*av_les| id.num + year
                                 , data = NoLow)
                 
                 # return the coefficients
                 coef(this_mod)[3]
               }
stopCluster(cl)

marg_les <- marginsplotdf2(model = les_mod, 
                          xterm = "treat2", 
                          zterm = "av_les", 
                          zseq = seq(min(NoLow$av_les, na.rm = T), 
                                     max(NoLow$av_les, na.rm = T), 
                                     by =.1))
moderator <- NoLow$av_les

les_plot <- inter_plot(df = marg_les, z = z, dydx = dydx,
                       x_breaks = c(0, 1, 2.5, 3),
                       x_labs = c("0 Legislative\nEffectiveness", 
                                  "1",
                                  "2.5 Legislative\nEffectiveness",
                                  " "),
                       x_axis_lab = "Legislative Effectiveness \n(Averaged Over Tenure)", 
                       y_axis_lab = "Marginal effect\nof Hiring MC",
                       plot_title = NULL,
                       x_coord = 2.5, y_coord = 0.2,
                        effect_lab = paste("Estimate =",
                                           round(coef(les_mod)[length(coef(les_mod))],2),
                                           "\n SE =",
                                           round(sd(res), 2)),
                        df2 = NoLow, moderator = "av_les",
                       ylims = c(-.5, 0.5),
                       xlims = c(-0.01,3))
ggdraw(les_plot)

ggsave(les_plot, 
       filename = here(dirname(this_dir),"replication", "images", "InterLES.eps"),
       device = cairo_ps,
       width = 5.5, height =5.5)




#--------------------------------------------------------------------------
# Appendix G Political Connections and Rule Changes

# Figure G.1: Limited Evidence of Rule Changes Driving the Effect


CTRL <- how(within = permute::Within(type = "series"), # time-series ordering ...
            #plots = Plots(strata = NoLow$id.num)
            blocks = NoLow$id.num) # ... within firm-clusters

#create 3000 firm-cluster-permuted treatment assignments

set.seed(185)
perm <- shuffleSet(NoLow$treat2, nset = 3000, control = CTRL)#each row of 'perm' represents a unique assignment


NoLow$return_on_assets <- NoLow$pretax_income/NoLow$DS_total_assets
NoLow$return_on_assets <- ifelse(is.infinite(NoLow$return_on_assets) == T, NA, 
                                 NoLow$return_on_assets)
NoLow$intensity <- I(NoLow$fixed_assets/NoLow$DS_total_assets)

NoLow$intensity <- ifelse(is.infinite(NoLow$intensity) == T, NA, 
                            NoLow$intensity)

inter_forms <- c("lead_tax ~ treat2*log(intensity+1 )| id.num + year",
                 "lead_tax ~ treat2*log(DS_total_assets+4)| id.num + year",
                 "lead_tax ~ treat2*log(fixed_assets+2057427)| id.num + year",
                 "lead_tax ~ treat2*return_on_assets| id.num + year",
                 "lead_tax ~ treat2*log(employees+1)| id.num + year")

ri_inter_forms <- c("lead_tax ~ new_treat*log(intensity+1 )| id.num + year",
                    "lead_tax ~ new_treat*log(DS_total_assets+4)| id.num + year",
                    "lead_tax ~ new_treat*log(fixed_assets+2057427)| id.num + year",
                    "lead_tax ~ new_treat*return_on_assets| id.num + year",
                    "lead_tax ~ new_treat*log(employees+1)| id.num + year")

sample_est <- rep(NA, length(inter_forms))
ri_sd <- rep(NA, length(inter_forms))

for(f in 1:length(inter_forms)){
  
  sample_est[f] <-coef( felm(as.formula(inter_forms[f]),
                             data = NoLow) )[3]
  
  
  
  cl <- makeCluster(detectCores()-1)
  registerDoParallel(cl) # register the cluster
  
  res <- foreach(i = 1:nrow(perm), 
                 .combine = "rbind", 
                 .packages = "lfe") %dopar% {
                   
                   crnt_perm <- perm[i,]
                   NoLow$new_treat <- NoLow$treat2[crnt_perm]
                   
                   this_mod <-felm(formula = as.formula(ri_inter_forms[f])
                                   , data = NoLow)
                   
                   # return the coefficients
                   coef(this_mod)[3]
                 }
  stopCluster(cl)
  
  ri_sd[f] <- sd(res)
  
}


#-----------------------
# spatial lag


used <- NoLow %>%
  drop_na(lead_tax, treat2, new_sec, id.num, year)


used$sector <- as.character(used$new_sec)
used$id.num <- as.character(used$id.num)
used$year <- as.character(used$year)
used$firm_year <- paste(used$id.num, used$year, sep = "_")

spat.sector <- t(combn(used$new_sec, m = 2))

w.id <- t(combn(used$id.num, m=2))
w.FY <- t(combn(used$firm_year, m=2))


spat.sector <- as.data.frame(spat.sector)
w.id <- as.data.frame(w.id)

names(spat.sector) <- c("sector_i", "sector_j")

weight <- ifelse(spat.sector[,1] == spat.sector[,2] &
                   w.FY[,1] != w.FY[,2], 1, 0)

w.id <- row.names(used)
w.id <- t(combn(w.id, m = 2))

w.edge <- data.frame(w.FY, as.numeric(as.character(weight)))
w.edge[,1] <- as.character(w.edge[,1])
w.edge[,2]<- as.character(w.edge[,2])
names(w.edge)[3]<- "weight"

w.mat <- as.matrix(w.edge)
w.mat[,1] <- as.character(w.mat[,1])
w.mat[,2]<- as.character(w.mat[,2])

w.mat2 <- graph.edgelist(w.mat[,1:2], directed = F)
E(w.mat2)$weight <- as.numeric(as.character(w.edge$weight))

panel.weight <- get.adjacency(w.mat2, attr='weight')
panel.weight <- as.matrix(panel.weight)
panel.weight <- mat2listw(panel.weight)


mat <- get.adjacency(w.mat2, attr='weight')
mat <- as.matrix(mat)


used$slx <- mat %*% used$treat2

slx_mod <- felm(lead_tax ~ slx | id.num + year | 0 | id.num, 
                data = used)


#identify clusters
CTRL <- how(within = Within(type = "series"), # time-series ordering ...
            blocks = used$id.num) # ... within firm-clusters
set.seed(185)
usedperm <- shuffleSet(used$treat2, nset = 3000, control = CTRL)#each row of 'perm' represents a unique assignment


cl <- makeCluster(detectCores()-1)
registerDoParallel(cl) # register the cluster

res <- foreach(i = 1:nrow(usedperm), 
               .combine = "rbind", 
               .packages = "lfe") %dopar% {
                 
                 
                 
                 crnt_perm <- usedperm[i,]
                 used$new_treat <- used$treat2[crnt_perm]
                 
                 used$slx_perm <- mat %*% used$new_treat
                 
                 
                 slx_mod <- felm(lead_tax ~ #treat2 + 
                                   slx_perm | id.num + year | 0 | id.num, 
                                 data = used)
                 
                 # return the coefficients
                 coef(slx_mod)[1]
               }
stopCluster(cl)

slx_res <- coef( felm(lead_tax ~ #treat2 + 
                        slx | id.num + year | 0 | id.num, 
                      data = used) )[1]
slx_sd <- sd(res)

sample_est <- rbind(c(sample_est, slx_res),
                    c(ri_sd, slx_sd))
sample_est <- t(sample_est)

sample_est <- as.data.frame(sample_est)
names(sample_est) <- c("pe", "se")

sample_est$spec <- c("intensity", "size", "fixed", "roa", "employees", "neighbor")
sample_est$order <- c(5:1,6)

sample_est$pe <- as.numeric(as.character(sample_est$pe))
sample_est$se <- as.numeric(as.character(sample_est$se))

alt <- ggplot(sample_est, aes(x = pe, y = reorder(spec,order))) +
  geom_point() +
  geom_errorbarh(aes(xmin = pe - (1.96*se), xmax = pe + (1.96*se)), 
                 height = 0, size = .25)+
  geom_errorbarh(aes(xmin = pe - (1.645*se), xmax = pe + (1.645*se)), 
                 height = 0, size = .75)+
  theme_classic() +
  geom_vline(xintercept = 0, lty=2) +
  #geom_vline(xintercept = -0.076942, lty=3) +
  scale_y_discrete(#breaks = 1:6,
    labels = c("Employees X \n Revolving Door    ",
               "Returns on Assets X \n Revolving Door     ",
               "Fixed Assets X \n Revolving Door    ",
               "Total Assets X \n Revolving Door    ",
               "Capital Intensity X \n Revolving Door     ",
               "Neighborhood Effect")) +
  labs(y = NULL,
       x = "Coefficient")

ggsave(alt, 
       filename = here(dirname(this_dir),"replication", "images", "AlternativeExplanations.eps"),
       device = cairo_ps,
       width = 5.5, height =6)

#-------------------------------------------------------------------------------------
# Appendix H Robustness to Staggered Treatment and Temporal Dynamics

# Figure H.2: Event Study Difference-in-Differences.

time.to.hire <- NoLow

time.to.hire$comp.id <-time.to.hire %>%
  group_by(id.num) %>%
  group_indices()

time.to.hire$year <- as.numeric(as.character(time.to.hire$year))

helper <- time.to.hire %>%
  dplyr::filter(treat2 == 1) %>%
  dplyr::group_by(id.num) %>%
  dplyr::mutate(firstocc = min(year)) %>%
  dplyr::filter(duplicated(id.num) == F) %>%
  dplyr::select(id.num, firstocc)

time.to.hire2 <- left_join(time.to.hire, helper, by = "id.num")

time.to.hire2$firstocc <- ifelse(is.na(time.to.hire2$firstocc) == T, 0, time.to.hire2$firstocc)


time.to.hire2 <- filter(time.to.hire2,
                        is.na(log_tax) == F &
                          is.na(year) == F &
                          is.na(comp.id) == F)

out <- att_gt(yname = "log_tax",
              gname = "firstocc",
              idname = "comp.id",
              tname = "year",
              xformla = ~ 1,  
              allow_unbalanced_panel = T,
              data = time.to.hire2,
              est_method = "dr",
              clustervars = "comp.id",
              bstrap = T, cband = T,alp=.1,
              print_details = T,
              control_group = "notyettreated",
              anticipation = 0
)

es <- aggte(out, type="dynamic", alp = 0.1,
            min_e = -6, # before -6 it drops below 100
            max_e = 4,  # choose four years like in PanelMatch analysis
            balance_e = 4) # balance on four treated years


tax_est <- data.frame(estimate = es$att.egt,
                      se = es$se.egt,
                      et = -6:4)

p_tax <- ggplot(tax_est, aes(x = et, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = estimate - 1.645*se,
                    ymax = estimate + 1.645*se), width = 0) +
  geom_hline(yintercept = 0, lty = 3)+
  theme_classic() +
  geom_vline(xintercept = -0.5, lty = 2) +
  scale_x_continuous(breaks = -6:4,
                     labels = c(-6:-1, "0\nMC Hired",
                                1:4)) +
  annotate(geom = "text", x = 2, y = 0.15,
           label = "Post-Treatment Period") +
  annotate(geom = "text", x = -3, y = 0.15,
           label = "Pre-Treatment Period") +
  labs(x = "Years Until/Since MC is Hired",
       y = "Difference-in-Difference Estimate",
       title = "A: Tax Rates")

p_tax


#-----------------------
# audit results

# prepare data for matching in pre-treatment outcome
helper <- time.to.hire2 %>%
  dplyr::group_by(comp.id) %>%
  dplyr::filter(year < firstocc) %>%
  dplyr::summarise(av_audit = mean(IRS_audit, na.rm = T)) %>%
  dplyr::mutate(av_audit = ifelse(is.nan(av_audit)==T, 0, av_audit))

# merge back in
time.to.hire3 <- left_join(time.to.hire2, helper, by = "comp.id")

time.to.hire3 <- as.data.frame(time.to.hire3)

time.to.hire3 <- time.to.hire3 %>%
  mutate(treated_firm = case_when(firstocc == 0 ~ 0,
                                  year < firstocc ~ 0,
                                  year >= firstocc ~ 1),
         lag_irs = dplyr::lag(IRS_audit, 1, order_by = year))

# remove missings, 
no_miss <- filter(time.to.hire3, year > 2006) 

out_aud <- att_gt(yname = "IRS_audit",
                  gname = "firstocc",
                  idname = "comp.id",
                  tname = "year",
                  xformla = ~ 1, 
                  panel = T,
                  allow_unbalanced_panel = T,
                  data = time.to.hire2,
                  est_method = "ipw",
                  clustervars = "comp.id",
                  bstrap = T, cband = T,alp=.1,
                  print_details = T,
                  control_group = "notyettreated",
                  anticipation = 0
)


#summary(out_aud)

#ggdid(out)

es_aud <- aggte(out_aud, type="dynamic", na.rm=T, alp = 0.1,
                min_e = -3, 
                max_e = 4,
                balance_e = 4)
#summary(es_aud)
#ggdid(es_aud)


audit_est1 <- data.frame(estimate = es_aud$att.egt,
                         se = es_aud$se.egt,
                         et = -3:4)

p_audit1 <- ggplot(audit_est1, aes(x = et, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = estimate - 1.645*se,
                    ymax = estimate + 1.645*se), width = 0) +
  geom_hline(yintercept = 0, lty = 3)+
  theme_classic() +
  geom_vline(xintercept = -0.5, lty = 2) +
  scale_x_continuous(breaks = -6:4,
                     labels = c(-6:-1, "0\nMC Hired",
                                1:4)) +
  annotate(geom = "text", x = 2, y = 0.6,
           label = "Post-Treatment Period") +
  annotate(geom = "text", x = -2, y = 0.6,
           label = "Pre-Treatment Period") +
  labs(x = "Years Until/Since MC is Hired",
       y = "Difference-in-Difference Estimate",
       title = "B: Audit Rates")

p_audit1

# remove small groups and missings
no_miss2 <- filter(time.to.hire3, is.na(lag_irs) == F & year > 2006)

out_aud2 <- att_gt(yname = "IRS_audit",
                   gname = "firstocc",
                   idname = "comp.id",
                   tname = "year",
                   xformla = ~ lag_irs,
                   panel = T,
                   allow_unbalanced_panel = T,
                   data = no_miss2,
                   est_method = "ipw",
                   clustervars = "comp.id",
                   bstrap = T, cband = T,alp=.1,
                   print_details = T,
                   control_group = "notyettreated",
                   anticipation = 0
)



es_aud2 <- aggte(out_aud2, type="dynamic", na.rm=T, alp = 0.1,
                 min_e = -3, 
                 max_e = 4,
                 balance_e = 4)

audit_est2 <- data.frame(estimate = es_aud2$att.egt,
                         se = es_aud2$se.egt,
                         et = -3:4)

p_audit2 <- ggplot(audit_est2, aes(x = et, y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin = estimate - 1.645*se,
                    ymax = estimate + 1.645*se), width = 0) +
  geom_hline(yintercept = 0, lty = 3)+
  theme_classic() +
  geom_vline(xintercept = -0.5, lty = 2) +
  scale_x_continuous(breaks = -6:4,
                     labels = c(-6:-1, "0\nMC Hired",
                                1:4)) +
  annotate(geom = "text", x = 1, y = 0.6,
           label = "Post-Treatment Period") +
  annotate(geom = "text", x = -2, y = 0.6,
           label = "Pre-Treatment Period") +
  labs(x = "Years Until/Since MC is Hired",
       y = "Difference-in-Difference Estimate",
       title = "C: Audit Rates (Matched)")

p_audit2


cs_plot <- plot_grid(p_tax, p_audit1, p_audit2)

ggsave(plot = cs_plot,
       filename = here(dirname(this_dir),"replication", "images", "cs_res.eps"),
       width = 11.1, height = 7,
       device = cairo_ps)

# Figure H.3: Robustness of Two-Way Fixed Effects.

# assume treatment for four years
NoLow$dummy_treat <- ifelse(NoLow$year ==  NoLow$primary.year.start.job |
                              NoLow$year ==  NoLow$primary.year.start.job +1 |
                              NoLow$year ==  NoLow$primary.year.start.job +2 | 
                              NoLow$year ==  NoLow$primary.year.start.job + 3, 1, 0)

NoLow$int.year <- as.integer(NoLow$year)
NoLow <- as.data.frame(NoLow)

NoLow$id.num <- as.integer(NoLow$id.num)


audit <- PanelMatch(lag = 1, 
                    time.id = "int.year", unit.id = "id.num", 
                    treatment = "dummy_treat", 
                    refinement.method = "mahalanobis", #alternative none
                    data = NoLow, 
                    match.missing = T, 
                    covs.formula = ~ I(lag(IRS_audit, 4)) + I(lag(cycle,2)), 
                    qoi = "att" ,
                    outcome.var = "lead_audit",
                    lead = 0, 
                    forbid.treatment.reversal = F)
audit_results <- PanelEstimate(sets = audit, 
                               data = NoLow, confidence.level = .9)
summary(audit_results)


bureau <- PanelMatch(lag =1, 
                     time.id = "int.year", unit.id = "id.num", 
                     treatment = "dummy_treat", 
                     refinement.method = "mahalanobis", #alternative none
                     data = NoLow, 
                     match.missing = T, 
                     covs.formula = ~ I(lag(log_tax, 4)) 
                     + I(lag(cycle,2)), 
                     qoi = "att" ,
                     outcome.var = "lead_tax",
                     lead = 0, 
                     forbid.treatment.reversal = T)
#bureau
bureau_obs <- length(print(bureau$att)[,3]) + sum(print(bureau$att)[,3])

bureau_results <- PanelEstimate(sets = bureau, 
                                data = NoLow, confidence.level = .9)

#summary(bureau_results)


pm_res <- rbind(data.frame(summary(audit_results), time = 1, DV = "Audits"),
                data.frame(summary(bureau_results), time = 1, DV = "Taxes"))
names(pm_res)[1:4] <- c("PE", "SE", "lwr90", "upr90")

p_dd <- ggplot(pm_res, aes(y = DV, x = PE)) +
  geom_point() +
  geom_errorbarh(aes(xmin = lwr90, xmax = upr90), height = 0) +
  geom_vline(xintercept = 0, lty = 3) +
  theme_classic() +
  labs(x = "Difference-in-Differences", y = NULL)

p_dd

ggsave(plot = p_dd,
       filename = here(dirname(this_dir),"replication", "images", "DD_robust.eps"),
       device = cairo_ps,
       width = 5, height = 4)


#--------------------------------------------------------------
# Appendix H.2 Event Study on Connected Boards

#Figure H.4: Event Study on Politically Connected Boards of Directors. 

treat_pan <- readRDS(file = here(dirname(this_dir), "replication", "data", "BX_evertreated_revolvers.rds"))

treat_pan$fyear <- as.integer(treat_pan$fyear)
treat_pan <- as.data.frame(treat_pan)
treat_pan <- subset(treat_pan, is.na(fyear) == F)

bureau <- PanelMatch(lag = 2, 
                     time.id = "fyear", unit.id = "anon_id", 
                     treatment = "treat", 
                     refinement.method = "mahalanobis", #alternative none
                     data = treat_pan, 
                     match.missing = T, 
                     covs.formula = ~ I(lag(etr, 4)) + I(lag(pi, 4)),  
                     qoi = "att" ,
                     outcome.var = "etr",
                     lead = 0:6, # try leads until sign-reversal
                     forbid.treatment.reversal = F, 
                     use.diagonal.variance.matrix = F)

set.seed(393939)
bureau_results <- PanelEstimate(sets = bureau, 
                                data = treat_pan, confidence.level = .9,
                                df.adjustment = T)

res <- as.data.frame(summary(bureau_results)$summary[1:7,c(1,3:4)])
res$time <- rownames(res)
res$order <- 1:7
names(res)[2:3] <- c("lwr","upr")

boots <- bureau_results$bootstrapped.estimates

mean_iter <- mean(apply(boots[,2:ncol(boots)], 1, mean))

boot_ci <- quantile(apply(boots[,2:ncol(boots)], 1, mean), probs = c(0.05, 0.95))

p1 <- ggplot(res, aes(x = reorder(time, order), y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin =  lwr,
                    ymax =  upr), width = 0) +
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  annotate(geom = "text", x = 4.5, y = 1,
           label = round(mean_iter, digits = 3))+
  geom_bracket(inherit.aes = FALSE,
               xmin = 1.75, xmax = 7.25, y.position = .6,
               label = paste0("[", round(boot_ci[1], digits = 3), "; ",
                              round(boot_ci[2], digits = 3), "]")
  ) +
  labs(x = "Time Since MC Hired",
       y = "Change in Effective Tax Rate\n(Difference-in-Differences Estimate)",
       title = "A: MC Arriving and Tax Rate") 

#----------
# leave
leave <- PanelMatch(lag = 2, 
                    time.id = "fyear", unit.id = "anon_id", 
                    treatment = "mc_leave", 
                    refinement.method = "mahalanobis", #alternative none
                    data = treat_pan, 
                    match.missing = T, 
                    covs.formula = ~ I(lag(etr, 4)) + I(lag(pi, 4)),  
                    qoi = "att" ,
                    outcome.var = "etr",
                    lead = 0:2, 
                    forbid.treatment.reversal = F)

set.seed(100)

leave_results <- PanelEstimate(sets = leave, 
                               data = treat_pan, confidence.level = .90,
                               df.adjustment = TRUE)


leave_res_g <- as.data.frame(summary(leave_results)$summary[1:3,c(1,3:4)])
leave_res_g$time <- rownames(leave_res_g)
leave_res_g$order <- 1:3
names(leave_res_g)[2:3] <- c("lwr","upr")


p2 <- ggplot(leave_res_g, aes(x = reorder(time, order), y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin =  lwr,
                    ymax =  upr), width = 0) +
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "Time Since MC Leaves Firm",
       y = "Change in Effective Tax Rate\n(Difference-in-Differences Estimate)",
       title = "B: MC Leaving and Tax Rate") 

#-----------------------
# audit data, but incomplete

treat_pan_audit <- readRDS(file = here(dirname(this_dir), "replication", "data", "audit_and_end_date.rds"))

treat_pan_audit <- as.data.frame(treat_pan_audit)
treat_pan_audit$fyear <- as.integer(treat_pan_audit$fyear)

audit <- PanelMatch(lag = 2, 
                    time.id = "fyear", unit.id = "anon_id", 
                    treatment = "treat", 
                    refinement.method = "mahalanobis", #alternative none
                    data = treat_pan_audit, 
                    match.missing = T, 
                    covs.formula = ~ I(lag(IRS_audit, 4)) + I(lag(pi, 4)), 
                    #size.match = 2, 
                    qoi = "att" ,
                    outcome.var = "IRS_audit",
                    lead = 0:3, 
                    forbid.treatment.reversal = F)

set.seed(8)

audit_results <- PanelEstimate(sets = audit, 
                               data = treat_pan_audit, confidence.level = .9,
                               df.adjustment = TRUE)


audit_res_g <- as.data.frame(summary(audit_results)$summary[1:4,c(1,3,4)])
audit_res_g$time <- rownames(audit_res_g)
audit_res_g$order <- 1:4
names(audit_res_g)[2:3] <- c("lwr", "upr")

p3 <- ggplot(audit_res_g, aes(x = reorder(time, order), y = estimate)) +
  geom_point() +
  geom_errorbar(aes(ymin =  lwr,
                    ymax =  upr), width = 0) +
  theme_classic() +
  geom_hline(yintercept = 0, lty = 3) +
  labs(x = "Time Since MC Hired",
       y = "Change in Audit Probability\n(Difference-in-Differences Estimate)",
       title = "C: Audit Probability") 


p <- plot_grid(p1, p2, p3)
p

ggsave(plot = p, 
       filename = here(dirname(this_dir),"replication", "images", "DD_BX_dates.eps"),
       height = 10, width = 8)


#-------------------------------------------------------------------
# Appendix I Heterogeneity by Position in Firm

# Figure I.5: Heterogeneity by the Revolver’s Position in the Firm.

NoLow$position <- tolower(NoLow$primary.position)

dir_index <- grep(paste(c("director", "board member", "chairman"),
                        collapse="|"), NoLow$position)

NoLow$director <- 0
NoLow$director[dir_index] <- 1


non_dir_board_index <- grep(paste(c("supervisory board", "board of advisors"),
                                  collapse="|"), NoLow$position)

NoLow$non_director <- 0
NoLow$non_director[non_dir_board_index] <- 1

lobby_index <- grep(paste(c("senior advisor", "consultant", "lobbyist",
                            "director of federal affairs"),
                          collapse="|"), NoLow$position)
NoLow$lobby <- 0
NoLow$lobby[lobby_index] <- 1

csuite_index <- grep(paste(c("vice-president", "ceo", "vice-president"),
                           collapse="|"), NoLow$position)
NoLow$csuite <- 0
NoLow$csuite[csuite_index] <- 1

tax_mod <- feols(lead_tax~  
                      treat2:director +
                      treat2:non_director + 
                      treat2:lobby +
                      treat2:csuite|id.num+year,
                    cluster = ~ id.num,
                    data = NoLow)

audit_mod <- feols(lead_audit ~  
                     treat2:director +
                     treat2:non_director + 
                     treat2:lobby +
                     treat2:csuite|id.num+year,
                   cluster = ~ id.num,
                   data = NoLow)

res <- rbind(data.frame(tidy(tax_mod), DV = "A: Tax Rate" ), 
             data.frame(tidy(audit_mod), DV = "B: Audit" ))

res$position <- c("Director", "Other Boards", "Gov't Affairs",
                  "C-Suite")

p<- ggplot(res, aes(x = estimate, y = position)) +
  geom_point() +
  geom_errorbarh(aes(xmin = estimate - 1.645*std.error,
                     xmax = estimate + 1.645*std.error),
                 height = 0)+
  theme_classic() +
  geom_vline(xintercept = 0, lty = 3) +
  facet_wrap(~ DV) +
  labs(x = "Estimate", y = NULL)

p

ggsave(plot = p,
       filename = here(dirname(this_dir),"replication", "images", "Het_by_position.eps"),
       device = cairo_ps,
       width = 8, height = 5)


#-------------------------------------------------------------------------
# Appendix J Revolving Door Contract Lobbyists

# Table J.1: Contracting with Revolver Lobbyists and Corporate Tax Rates


rev_tax_mod <- felm(lead_tax ~ revolver | id.num + year|0|id.num,
                    data = NoLow)

rev_audit_mod <- felm(lead_audit ~ revolver | id.num + year|0|id.num,
                      data = NoLow)


stargazer(rev_tax_mod, rev_audit_mod, 
          covariate.labels = "Contract w. Revolver", omit.stat = c("rsq","adj.rsq",  "f"), 
          df =F, dep.var.labels = c("ln Tax Rate", "IRS Audit"),
          add.lines = list(c("Firm FE?", "Yes", "Yes", "Yes", "Yes", "Yes"),
                           c("Year FE?", "Yes", "Yes", "Yes", "Yes", "Yes")),
          title = c("Contracting with Revolver Lobbyists and Corporate Tax Rates"),
          label = c("tab:contract"),
          out = here(dirname(this_dir),"replication", "tables", "contract_with_revolver.tex"))


#.-------------------------------------------------------------
# Appendix K Sensitivity Analysis

# Figure K.6: Tipping Point Analysis of Effective Tax Rate. 


mod <- felm(lead_tax ~ treat2|id.num+year|0|0,
            data = NoLow)


audit_mod <- felm(lead_audit ~ treat2|id.num+year|0|0,
                  data = NoLow)


#--------------------------
# baseline estimates

# other strategies: lobbying
helper <- felm(log_expend ~ treat2|id.num+year|0|id.num,
               data = NoLow)
summary(helper)

out_ass <- felm(lead_tax ~ log_expend|id.num+year|0|0,
                data = NoLow)

coef_lobby <- coef(out_ass)[1]

#

lobby_smd <- tidy(helper) %>% filter(term == "treat2") %>% pull("estimate") /glance(helper) %>% pull("sigma")

# other strategies: donations
helper_don <- felm(donate ~ treat2|id.num+year|0|id.num,
                   data = NoLow)

out_ass_don <- felm(lead_tax ~ donate|id.num+year|0|0,
                    data = NoLow)

coef_don <- coef(out_ass_don)[1]

don_smd <- tidy(helper_don) %>% filter(term == "treat2") %>% pull("estimate") /glance(helper) %>% pull("sigma")


# firm differences: finances

covar <- c("lag_capital",
           "lag_assets", 
           "lag_ent",
           "lag_emp", 
           "lag_sales", 
           "lag_income", 
           "lag_turnover",   
           "lag_market_value",   
           "lag_prices")

forms <- paste0(covar, "~ treat2|id.num + year")
forms2 <- paste0("lead_tax ~", covar, "|id.num+year")
forms3 <- paste0("lead_audit ~", covar, "|id.num+year")

comp_est <- list()

audit_est <- list()

for(i in 1:length(forms)){
  
  crnt_mod <- felm(formula = as.formula(forms[i]),
                   data = NoLow)
  est_sigma <- glance(crnt_mod) %>% pull("sigma")
  
  out_mod <- felm(formula = as.formula(forms2[i]),
                  data = NoLow)
  
  audit_out_mod <- felm(formula = as.formula(forms3[i]),
                        data = NoLow)
  audit_sigma <- glance(audit_out_mod) %>% pull("sigma")
  
  comp_est[[i]] <- data.frame(smd = coef(crnt_mod)/est_sigma,
                              out = coef(out_mod),
                              covar = covar[i]) 
  
  audit_est[[i]] <- data.frame(smd = coef(crnt_mod)/est_sigma,
                               out = coef(audit_out_mod),
                               std_out = coef(audit_out_mod)/audit_sigma,
                               covar = covar[i]) 
}

comp_est <- do.call("rbind", comp_est)
audit_est <- do.call("rbind", audit_est)


#--------------------------
# get ready for sensitivity
res <- tidy(mod) %>% filter(term == "treat2")
res$std.error <- corr_se
res$conf.low = res$estimate - 1.645*res$std.error
res$conf.high = res$estimate + 1.645*res$std.error

res <- res %>% select(-c(statistic, p.value))

plot_stuff <- res %>%
  select(conf.high) %>%
  tip_coef(exposure_confounder_effect = seq(0.01, .25, by = 0.05),
      verbose = TRUE)

p1<-ggplot(plot_stuff, aes(x = exposure_confounder_effect, 
                           y = confounder_outcome_effect)) +
  geom_line(color = "black") +
  scale_y_continuous(limits = c(-2, 0.1)) +
  labs(x = "Standardized mean difference\nbtwn Hiring and not-yet-hiring firms",
       y = "Association with ln ETR",
       title = "A: Single Confounders Leading to Insignificance") +
  theme_bw() +
  annotate(geom = "text", x = lobby_smd, y = coef_lobby+0.1,
           label = "Lobbying") +
  geom_point(inherit.aes = FALSE,
             aes(x = lobby_smd, y = coef_lobby)) +
  annotate(geom = "text", 
           x = comp_est$smd[grep("lag_asset", comp_est$covar)], 
           y = comp_est$out[grep("lag_asset", comp_est$covar)]+0.1,
           label = "Total Assets") +
  geom_point(inherit.aes = FALSE,
             aes(x = comp_est$smd[grep("lag_asset", comp_est$covar)], 
                 y = comp_est$out[grep("lag_asset", comp_est$covar)]))+
  annotate(geom = "text", 
           x = comp_est$smd[grep("lag_income", comp_est$covar)]-0.024, 
           y = comp_est$out[grep("lag_income", comp_est$covar)]+0.1,
           label = "Gross Income") +
  geom_point(inherit.aes = FALSE,
             aes(x = comp_est$smd[grep("lag_income", comp_est$covar)], 
                 y = comp_est$out[grep("lag_income", comp_est$covar)]))

p1

# number of confounders necessary

# for many values
confound <- seq(0.03, 0.2, by = 0.01)

n_con <- list()

for(i in 1:18){
  
  n_con[[i]] <- res %>%
    select(conf.high) %>%
    tip_coef(exposure_confounder_effect = confound[i],
             confounder_outcome_effect = -confound[i],
             verbose = TRUE)

}

n_con <- do.call("rbind", n_con)

p2<-ggplot(n_con, aes(x = exposure_confounder_effect,
                      y = n_unmeasured_confounders)) +
  geom_line() +
  theme_bw()+
  annotate(geom = "text", x = lobby_smd+0.005, y = 4,
           label = "Lobbying") +
  geom_point(inherit.aes = FALSE,
             aes(x = lobby_smd, y = 3)) +
  labs(x = "Standardized mean difference and\nassociation with ETR",
       y = "Number of confounders needed\ntorender estimate insignificant",
       title = "B: Number of Similar Confounders")


p <- plot_grid(p1, p2)
p

ggsave(plot = p,
       filename = here(dirname(this_dir),"replication", "images", "sensitivity.eps"),
       device = cairo_ps,
       width = 11,
       height = 6)

# Figure K.7: Tipping Point Analysis of Audit Initiation.
audit_res <- tidy(mod2) %>% filter(term == "treat2")
audit_res$std.error <- corr_se_audit
audit_res$conf.low = audit_res$estimate - 1.645*audit_res$std.error
audit_res$conf.high = audit_res$estimate + 1.645*audit_res$std.error

audit_res <- audit_res %>% select(-c(statistic, p.value))

audit_plot_stuff <- audit_res %>%
  select(conf.high) %>%
  tip_coef(exposure_confounder_effect = seq(0.1, .5, by = 0.05),
           verbose = TRUE)

# correlation lobby and audit

audit_out_ass <- felm(lead_audit ~ lobby|id.num+year|0|0,
                      data = NoLow)
coef_audit_lobby <- coef(audit_out_ass)[1]


p3<-ggplot(audit_plot_stuff, aes(x = exposure_confounder_effect,
                                 y = confounder_outcome_effect)) +
  geom_line(color = "black") +
  scale_y_continuous(limits = c(-1.05, 0.3)) +
  labs(x = "Standardized mean difference\nbtwn Hiring and not-yet-hiring firms",
       y = "Association with Audit Initiation",
       title = "A: Single Confounders Leading to Insignificance") +
  theme_bw() +
  annotate(geom = "text", 
           x = audit_est$smd[grep("lag_asset", audit_est$covar)], 
           y = audit_est$out[grep("lag_asset", audit_est$covar)]+0.1,
           label = "Total Assets") +
  geom_point(inherit.aes = FALSE,
             aes(x = audit_est$smd[grep("lag_asset", audit_est$covar)], 
                 y = audit_est$out[grep("lag_asset", audit_est$covar)]))+
  annotate(geom = "text", 
           x = audit_est$smd[grep("lag_ent", audit_est$covar)]+0.06, 
           y = audit_est$std_out[grep("lag_ent", audit_est$covar)]+0.15,
           label = "Enterprise Value") +
  geom_point(inherit.aes = FALSE,
             aes(x = audit_est$smd[grep("lag_ent", audit_est$covar)], 
                 y = audit_est$std_out[grep("lag_ent", audit_est$covar)]))

p3


confound <- seq(0.01, 0.5, by = 0.02)

n_con <- list()

for(i in 1:length(confound)){
  
  n_con[[i]] <-audit_res %>%
    select(conf.high) %>%
    tip_coef(exposure_confounder_effect = confound[i],
             confounder_outcome_effect = -0.1,
             verbose = TRUE)
}

n_con <- do.call("rbind", n_con)

p4<-ggplot(n_con, aes(x = exposure_confounder_effect, y = n_unmeasured_confounders)) +
  geom_line() +
  theme_bw()+
  labs(x = "Standardized mean difference",
       y = "Number of confounders needed\ntorender estimate insignificant",
       title = "B: Number of Similar Confounders")

p4

p5 <- plot_grid(p3, p4)
p5


ggsave(plot = p5,
       filename = here(dirname(this_dir),"replication", "images", "audit_sensitivity.eps"),
       device = cairo_ps,
       width = 11,
       height = 5)


